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Abstract 

The  least  absolute  deviations  criterion,  or  the  t\  norm,  is  frequently  used  for  approx- 
imation where  the  data  may  contain  outliers  or  ‘wild  points’.  One  of  the  most  popular 
methods  for  solving  the  least  absolute  deviations  data  fitting  problem  is  the  Barrodale 
and  Roberts  (BR)  algorithm  (1973),  which  is  based  on  linear  programming  techniques 
and  the  use  of  a modified  simplex  method  [1].  This  algorithm  is  particularly  efficient. 
However,  since  it  is  based  upon  the  simplex  method  it  can  be  susceptible  to  the  accu- 
mulation of  unrecoverable  rounding  errors  caused  by  using  an  inappropriate  pivot.  In 
this  paper  we  shall  show  how  we  can  extend  a numerically  stable  form  of  the  simplex 
method  to  the  special  case  of  i\  approximation  whilst  still  maintaining  the  efficiency  of 
the  Barrodale  and  Roberts  algorithm.  This  extension  is  achieved  by  using  the  i\  char- 
acterization to  rebuild  the  relevant  parts  of  the  simplex  tableau  at  each  iteration.  The 
advantage  of  this  approach  is  demonstrated  most  effectively  when  the  observation  matrix 
of  the  approximation  problem  is  sparse,  as  in  the  case  when  using  compactly  supported 
basis  functions  such  as  B-splines.  Under  these  circumstances  the  new  method  is  consid- 
erably more  efficient  than  the  Barrodale  and  Roberts  algorithm  as  well  as  being  more 
robust. 


1 Introduction 


Given  a set  of  m data  points  { (x, , t/j)  }£Lj , the  , or  least  absolute  deviations  curve-fitting 
problem  seeks  c € ]Rn  to  solve  the  optimization  problem 


min  \\y-Ac\\i 

C 


E 


n 

Vi  ^ ^ a,,  j Cj 
3= 1 


(1.1) 


where  A is  an  m x n observation  matrix,  and  r,  denotes  the  residual  of  the  ith  point. 

Another  way  of  stating  the  i\,  or  least  absolute  deviations  curve-fitting  problem,  is  by 
the  characterization  theory  of  an  ii  solution  [8],  which  may  be  given  in  different  forms. 
The  following  is  perhaps  the  most  commonly  used. 
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A vector  c 6 IRn  solves  the  minimization  problem  (1.1)  if  and  only  if  there  exist 
A € Itm  such  that 


AtA  - 0 with 


f |A*|  < 1,  for  i e Z, 
\ A i = sign(j-j),  for  i £ Z , 


(1.2) 


where  Z represents  the  set  of  indices  for  which  r,  = 0. 

One  of  the  popular  methods  designed  for  solving  the  l\  approximation  problem  is  the 
Barrodale  and  Roberts  (BR)  algorithm.  It  replaces  the  unconstrained  variables  c and  r 
in  (1.1)  by  nonnegative  variables  c+,  c~,  u and  v,  and  considers  the  linear  programming 
problem 


min  eTu  + eTv 

C 

subject  to  Ac+  — Ac~  + u — v = y,  (1.3) 

c+,c~,u,v  > 0. 


Much  of  the  reason  for  the  popularity  of  the  BR  algorithm  is  that  it  exploits  the 
characteristics  of  the  t\  approximation  in  order  to  solve  the  problem  in  a more  efficient 
manner  than  the  general  simplex  approach.  However,  it  is  a simplex  based  method,  and 
so  it  is  susceptible  to  numerical  instabilities  caused  by  using  inappropriate  pivots.  The 
new  method  presented  here  uses  matrix  factorization  instead  of  simplex  pivoting.  This 
approach  allows  numerically  stable  updates  to  be  made,  thus  avoiding  the  unnecessary 
build-up  of  rounding  errors.  This  method  is  particularly  efficient  when  the  observation 
matrix  is  large  and  sparse  [5]. 

Bartels  [2]  and  Gill  and  Murray  [4]  presented  methods  that  concentrate  on  avoiding 
the  inherent  instability  of  the  simplex  method.  However,  these  methods  are  designed  for 
a general  linear  programming  problem  and  if  we  were  to  employ  these  techniques  for  the 
special  case  of  the  l\  problem,  the  storage  requirements  and  computational  workload  of 
the  method  would  be  unnecessarily  large  compared  to  those  of  the  highly  efficient  BR 
algorithm. 

The  l\  problem  is,  in  essence,  an  interpolation  problem.  The  aim  of  any  iterative 
procedure  for  the  t\  problem  is  to  find  an  optimal  set  of  interpolation  points.  Indeed, 
this  is  how  the  BR  algorithm  solves  the  i\  problem.  It  begins  with  all  coefficients,  c, 
set  to  zero  (being  non-basic  variables),  and  during  each  iteration  of  stage  one,  one  of 
the  residuals,  r-j,  becomes  non-basic  by  making  the  corresponding  point  an  interpolation 
point  (i.e.,  the  coefficients  are  altered  so  that  Vi  = 0).  At  the  end  of  stage  one,  the  current 
estimate  interpolates  n distinct  points.  During  stage  two,  the  interpolation  points  are 
exchanged  one  at  a time  with  a non-interpolation  point  until  an  optimal  solution  is 
achieved. 

In  fact,  the  new  algorithm  is  effectively  identical  to  the  BR  algorithm  in  the  sense  that 
we  use  exactly  the  same  pivoting  strategy.  However,  we  start  with  a predetermined  set 
of  interpolation  points  and  do  not  store  the  simplex  tableau  directly.  In  each  iteration, 
we  only  reconstruct  the  parts  of  the  simplex  tableau  that  are  needed  by  the  more  stable 
approach  employed. 
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2 A more  stable  computational  approach 

The  linear  programming  presentation  of  a least  absolute  deviations  curve-fitting  problem 
is  given  in  (1.3).  It  is  a standard  linear  programming  problem  of  dimension  mx  (2m+2n). 
The  robust  approaches  of  Bartels  and  Gill  and  Murray  can  be  applied  to  solve  it.  They 
involve  the  factorization  of  an  m x m matrix.  On  the  other  hand,  the  BR  algorithm  only 
deals  with  an  m x n matrix  in  each  iteration,  if  m S>  n,  the  direct  usage  of  these  stable 
approaches  is  less  efficient.  We  shall  show  next  that  the  factorization  of  an  n x n matrix 
is  all  that  is  required  at  each  iteration. 

We  split  the  data  points  based  on  the  set  interpolation  Z,  and  let  Az,  yz , and 
vz  be  the  counterparts  of  A,  y,  u and  v in  (1.3)  corresponding  to  the  set  Z.  Their 
complementary  matrix  and  vectors  are  denoted  by  Az  and  yz,  uz  and  vz,  so  that  Az 
and  Az  comprise  A,  etc.,  problem  (1.3)  can  be  expressed  as 

min  eT(uz + uz)  + eT(vz + vz) 

C 

subject  to  Azc+  — Azc~  + uz  - vz 
Azc+  - Azc~  +uz  - vz 
c+,c~,uz,uz,vz,vz 

Since  the  coefficients  for  cj  are  just  the  negative  of  the  coefficients  for  Cj  , j = 1,2 
it  is  possible  to  suppress  cj  and  let  c represent  the  unconstrained  variable.  The  initial 
simplex  tableau  associated  with  problem  (2.1)  can  be  constructed  in  matrix  form  by 
Table  1,  where  e*,,  k = m,n,m  - n,  are  1x1  vectors  with  all  components  equal  to  one. 


BV 

c 

Uz 

uz 

vz 

vz 

r 

uz 

I 

0 

-I 

0 

yz 

Uz 

Az 

0 

I 

0 

-I 

yz 

Z 

*(&) 

0 

0 

-2  el 

-2eT 

^m-n 

6t  ( 1*  ^ 
m\yz  J 

Tab.  1.  The  initial  simplex  tableau  of  the  i\  fitting  problem. 


As  we  know,  the  simplex  method  is  an  iterative  procedure  in  which  each  iteration  is 
characterized  by  specifying  which  m of  2m +n  variables  are  basic.  For  the  l\  approxima- 
tion, we  are  only  concerned  with  those  vertices  which  are  formed  by  a set  of  interpolation 
points.  For  n interpolation  points,  the  basic  variables  consist  of  n of  the  coefficient  para- 
meters c and  m — n of  the  parameters  uz  corresponding  to  the  non-interpolation  points. 

Let  B be  the  m x to  basis  matrix  whose  columns  consist  of  the  m columns  associated 
with  the  basic  variables.  Then 


= yz, 
= Vz, 

> o. 


(2.1) 


A robust  algorithm,  for  least  absolute  deviations  curve  fitting  473 


BV 

uz 

r 

c 

A-1 

Az 

A~zlyz 

Uz 

-AzA~i 

rz 

Z 

-el_n(AzA-z')-el 

^m—nTZ 

Tab.  2.  The  condensed  simplex  tableau  associated  with  a set  of  interpolation 
points. 


It  is  readily  verified  that  the  inverse  of  B can  be  written  in  the  form  of  (2.3)  as  long 
as  Az  is  invertible. 


Equation(2.3)  shows  that  the  explicit  inverse  computation  of  an  m x m matrix  in  the 
form  of  (2.2)  can  be  achieved  by  dealing  with  an  inverse  of  an  n x n matrix,  and  in 
general,  n <C  m. 

To  make  the  m non-basic  variables  become  basic,  we  multiply  the  whole  simplex 
tableau  by  B~l,  and  omit  the  identity  and  zero  matrices.  Then  new  simplex  tableau  is 
given  in  Table  2. 

An  arbitrary  choice  of  the  interpolation  set  Z may  cause  some  of  the  values  in  the 
right  hand  side  column  to  become  negative.  Although  it  is  permissible  for  the  coefficient 
parameters  c to  be  negative,  for  those  rows  having  negative  residuals  rz,  we  restore 
feasibility  by  exchanging  the  corresponding  uz  for  vz-  This  exchanging  can  be  made  by 
subtracting  twice  those  rows  from  the  objective  row  and  changing  the  sign  of  the  original 
rows  [1]. 

Such  an  exchange  process  can  be  expressed  in  matrix  terms  by  introducing  a sign 
vector 

A z = sign  (rz). 

Let  Azs  represent  the  matrix  which  is  obtained  by  multiplying  those  rows  of  Az  asso- 
ciated with  negative  residuals  by  —1, 


Azs  = diag(Az)Az- 
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BV 

Uz 

r 

c 

Az1 

A~zlyz 

Uz 

-Azs*z 

Vz\ 

Z 

-X  l(AzA~zl)-el 

X\rz 

Tab.  3.  Restoration  of  feasibility  of  the  simplex  tableau. 

The  simplex  tableau  after  restoring  feasibility  is  shown  in  Table  3. 

The  point  to  be  removed  from  Z is  decided  by  the  values  of  the  objective  row.  Each 
time  the  maximum  value  of  the  objective  row  (including  the  suppressed  columns)  is 
chosen,  we  let  the  index  of  this  element  be  k.  In  order  to  choose  which  new  point  is 
to  join  the  set  Z,  we  compute  the  value  of  the  pivotal  column,  the  fcth  column  in  the 
simplex  tableau.  Since  the  simplex  tableau  is  in  the  form  of 

A-1 

the  fcth  column  can  be  obtained  by  using  Azs  and  the  fcth  column  of  A],1 . 

The  BR  algorithm  pivoting  strategy  is  adopted  to  decide  which  new  point  is  to  be 
added  to  the  interpolation  set,  when  a new  set  of  indices  Z is  generated.  We  repeat  the 
process  in  an  iterative  manner  until  the  optimal  solution  is  achieved. 

Table  3 is  in  fact  identical  to  the  simplex  tableau  of  the  BR  algorithm  in  stage  2. 
The  difference  here  is  that  the  BR  algorithm  is  implemented  by  a simplex  pivoting 
approach,  while  the  transformation  of  the  simplex  tableau  in  the  form  of  Table  3 can  be 
accomplished  in  a numerically  more  stable  manner. 

3 The  improved  method 

The  improved  method  starts  with  a predetermined  interpolation  set  Z,  the  minimum 
requirement  for  Z being  that  it  forms  a well-behaved  matrix  Az-  For  B-spline  basis  func- 
tions, we  can  choose  any  set  of  points  satisfying  the  Schoenberg- Whitney  condition  [6] . 
For  a Chebyshev  polynomial  basis,  points  close  to  the  n Chebyshev  zeros  can  be  regarded 
as  the  initial  interpolation  set.  In  other  cases,  we  can  choose  points  approximate  to  them 
or  even  uniformly  distributed. 

If  we  denote  the  set  of  A*,  i € Z,  as  \z,  we  can  rewrite  the  characterization  equa- 
tion (1.2)  as 

A"zXz  — —A'gXz , (3-1) 

and  Xz  can  be  obtained  mathematically  from 
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A z = -(^I)_1(^I^z)- 

(3.2) 

Table  3 shows  that  the  objective  row  can  be  computed  as 

Objective  row  = —(XzAz)A^1  — ejt. 

(3.3) 

Thus,  using  (3.2)  we  conclude  that 

Objective  row  = A^  - e„. 

(3.4) 

We  know  that  at  the  solution  all  the  values  in  the  objective  row  are  in  the  range 
[-2,  0],  and  also  |A|  < 1.  This  latter  result  can  be  explained  in  terms  of  the  former  by 
the  relationship  (3.4). 

(3.4)  is  useful  because  it  can  be  used  to  verify  whether  an  interpolation  set  forms  an 
optimal  solution,  or  to  compute  A from  the  values  of  the  objective  row.  We  use  it  to 
compute  the  values  of  the  objective  row. 

The  improved  method  can  be  summarized  as  follows; 

(1)  Choose  an  initial  set  of  interpolation  points  and  form  the  set  Z. 

(2)  Construct  Az,  yz  and  their  counterpart  Az,  yz  accordingly. 

(3)  Solve  the  equation  Azc  — yz  for  c,  and  compute 

rz  = yz~  Azc,  and  Xz  = sign(rz). 

(4)  Obtain  the  values  of  A z from  the  equation 

A^Xz  — —A'zXz-  (3-5) 

(5)  If  \Xz\  < 1 hold,  the  current  solution  is  optimal,  and  the  algorithm  terminates. 
Otherwise,  continue. 

(6)  Obtain  the  objective  row  of  the  BR  simplex  tableau  from 

objective  row  = A^  — e„. 

(7)  Examine  the  values  of  the  objective  row;  the  point  associated  with  the  maximum 
value  of  the  objective  row  is  chosen  to  leave  the  set  Z. 

(8)  Decide  the  point  to  add  by  the  BR  pivoting  strategy.  Obtain  a new  set  of  indices 
Z , and  repeat  from  step  2. 

4 Practical  considerations  and  application  to  the  t\  spline  ap- 
proximation 

The  robustness  of  the  above  algorithm  stems  from  the  reliable  updating  of  the  relevant 
parts  of  the  simplex  tableau  in  each  iteration.  The  major  computational  work  is  obtaining 
(explicitly  or  implicitly)  the  inverse  of  an  n x n matrix  AZ-  It  can  be  calculated  and 
stored  explicitly  by  using  an  LU  or  QR  factorization,  or  preferably  it  can  be  expressed  as 
a product  of  factors.  Since  Az  differs  from  its  predecessor  by  only  one  row,  savings  can 
be  made  by  reusing  results  from  the  previous  step.  Necessary  material  is  available  [4,  7] 
regarding  the  stable  implementation  of  this  row  updating  procedure. 
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m = 512 

Numbers  of  iterations 

Execution  Time  (seconds) 

Q = 

New 

BR 

New 

BR 

44 

57 

125 

1.6 

14.7 

49 

75 

111 

2.2 

13.4 

54 

71 

134 

2.4 

20.2 

59 

83 

156 

3.0 

26.8 

64 

78 

160 

3.1 

32. ‘ 

69 

88 

194 

4.0 

42.4 

74 

75 

165 

3.7 

36.0 

79 

87 

189 

4.8 

48.1 

Tab.  4.  The  number  of  iterations  and  execution  time  taken  by  the  algorithm  of 
this  paper  and  the  Barrodale  and  Roberts  algorithm  for  a set  of  512  response  data 
points  provided  by  the  National  Physical  Laboratory. 


Sparsity  almost  always  is  more  important  than  matrix  dimension.  Additional  savings 
can  be  made  if  the  observation  matrix  A is  sparse  or  structured.  Approximation  using 
a B-spline  basis  often  occurs  in  practical  applications.  In  such  cases,  A is  block  banded, 
and  Az  can  be  triangularized  using  0(n)  flops  [3],  Similarly,  the  sparsity  of  A can  be 
exploited  to  compute  other  relevant  parts  of  the  simplex  tableau  efficiently. 

We  have  applied  our  method  to  solve  the  least  absolute  deviations  curve-fitting  prob- 
lems by  B-splines  using  various  numbers  of  interior  knots.  All  software  was  written  in 
MATLAB  and  implemented  on  a Sun  Workstation.  The  initial  interpolation  points  are 
chosen  to  be  those  points  corresponding  to  the  maximum  value  in  each  column  of  the 
observation  matrix  A. 

Some  of  our  computational  results  are  reported  in  Tables  4 and  5.  Each  table  presents 
the  outcomes  of  a particular  set  of  data  points  by  the  new  method  and  by  the  BR 
algorithm. 

All  the  experimental  results  exhibit  the  effectiveness  of  the  improved  method  on  large, 
sparse  systems.  Although  these  tables  show  that  the  improved  method  is  faster  than  the 
BR  algorithm,  it  would  be  unfair  to  judge  the  convergence  speed  purely  based  upon 
the  time  taken,  since  the  improved  method  embodies  some  MATLAB  built-in  functions, 
while  the  BR  algorithm  uses  only  user-defined  functions.  However,  on  average,  the  new 
method  requires  far  fewer  iterations  than  the  BR  algorithm,  and  is  competitive  with  the 
BR  algorithm  both  in  efficiency  and  accuracy  for  a structured  system. 

Further  work  to  be  addressed  by  the  authors  will  involve  a definitive  implementation 
of  this  algorithm  in  Fortran,  and  development  of  an  error  analysis  for  both  the  improved 
method  and  the  BR  algorithm. 
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m = 1200 

Numbers  of  iterations 

Execution  Time  (seconds) 

9 = 

New 

BR 

New 

BR 

50 

82 

143 

4.0 

58.7 

56 

105 

165 

5.2 

85.8 

62 

113 

190 

6.1 

110.2 

68 

131 

189 

7.6 

110.4 

74 

121 

223 

7.8 

157.9 

80 

132 

216 

9.2 

163.2 

86 

155 

245 

11.8 

209.8 

92 

173 

252 

14.0 

241.8 

98 

153 

272 

13.6 

292.6 

Tab.  5.  The  number  of  iterations  and  execution  time  taken  by  the  algorithm 
of  this  paper  and  the  Barrodale  and  Roberts  algorithm  for  a set  of  1200  data 
points,  generated  by  MATLAB  command  x = linspace(l,  10, 1200)';  y = log(x)  + 
randn(1200, 1). 
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